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Abstract 

This article addresses some issues related to statistical description of gravitating 
systems in an expanding backgrounds. In particular, I describe (a) how the non 
linear mode-mode coupling transfers power from one scale to another in the Fourier 
space if the initial power spectrum is sharply peaked at a given scale and (b) what 
are the asymptotic characteristics of gravitational clustering that are independent 
of the initial conditions. The analysis uses a new approach based on an integro- 
differential equation for the evolution of the gravitational potential in the Fourier 
space. I show how this equation allows one to understand several aspects of nonlinear 
gravitational clustering and provides insight in to the transfer of power from one 
scale to another through nonlinear mode coupling. Numerical simulations as well 
as analytic work shows that power transfer leads to a universal power spectrum at 
late times, somewhat reminiscent of the existence of Kolmogorov spectrum in fluid 
turbulence. 
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1 Introduction 


One of the central problems in cosmology is to describe the non linear phases 
of the gravitational clustering starting from an initial spectrum of density 
fluctuations. It is often enough (and necessary) to use a statistical description 
and relate different statistical indicators (like the power spectra, nth order 
correlation functions etc.) of the resulting density distribution to the statisti¬ 
cal parameters (usually the power spectrum) of the initial distribution. The 
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relevant scales at which gravitational clustering is non linear are less than 
about 10 Mpc (where 1 Mpc = 3 x 10 24 cm is the typical separation between 
galaxies in the universe) while the expansion of the universe has a characteris¬ 
tic scale of about few thousand Mpc. Hence non linear gravitational clustering 
in an expanding universe can be adequately described by Newtonian gravity 
provided the rescaling of lengths due to the background expansion is taken 
into account. 


As to be expected, cosmological expansion introduces several new factors into 
the problem when compared with the study of statistical mechanics of isolated 
gravitating systems. (For a general review of statistical mechanics of gravitat¬ 
ing systems, see JI]. For a sample of different approaches, see [§] and the 
references cited therein. Review of gravitational clustering in expanding back¬ 
ground is also available in several textbooks in cosmology jl, 4].) Though this 
problem can be tackled in a ‘practical’ manner using high resolution numerical 
simulations (for a review, see [fj]), such an approach hides the physical princi¬ 
ples which govern the behaviour of the system. To understand the physics, it 
is necessary to attack the problem from several directions using analytic and 
semi analytic methods. Several such attempts exist in the literature based on 
Zcldovich(like) approximations 16], path integral and perturbative techniques 
[7], nonlinear scaling relations [§] and many others. In spite of all these it is 
probably fair to say that we still do not have a clear analytic grasp of this 
problem, mainly because each of these approximations have different domains 
of validity and do not arise from a central paradigm. 

I propose to attack the problem from a different angle, which has not received 
much attention in the past. The approach begins from the dynamical equation 
for the the density contrast in the Fourier space and casts it as an integro- 
differential equation. Though this equation is known in the literature (see, 
e.g. 0), it has received very little attention because it is not ‘closed’ mathe¬ 
matically; that is, it involves variables which are not natural to the formalism 
and thus further progress is difficult. I will, however, argue that there exists a 
natural closure condition for this equation based on Zcldovich approximation 
thereby allowing us to write down a closed integro-differential equation for the 
gravitational potential in the Fourier space. 


It turns out that this equation can form the basis for several further investi¬ 
gations some of which are described in ref. |s] and in the second reference in 
|l]. Here I will concentrate on just two specific features, centered around the 
following issues: 


• If the initial power spectrum is sharply peaked in a narrow band of wave¬ 
lengths, how does the evolution transfer the power to other scales? In par¬ 
ticular, does the non linear evolution in the case of gravitational interactions 
lead to a universal power spectrum (like the Kolmogorov spectrum in fluid 


2 


turbulence)? 

• What is the nature of the time evolution at late stages? Does the gravita¬ 
tional clustering at late stages wipe out the memory of initial conditions 
and evolve in a universal manner? 


2 An integral equation to describe nonlinear gravitational cluster¬ 
ing 


The expansion of the universe sets a natural length scale (called the Hubble 
radius) du = c(a/a )“ 1 which is about 4000 Mpc in the current universe. In 
any region which is small compared to du one can set up an unambiguous 
coordinate system in which the proper coordinate of a particle r (t) = a(t)x(t) 
satisfies the Newtonian equation r = —Vr® where <f> is the gravitational 
potential. The Lagrangian for such a system of particles is given by 
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we note that: (i) the total time derivative can be ignored; (ii) using a = 
— ( 47 rG/ 3 )pf,a, the term <L frw = —(l/2)aaxf = (2ttG / 3) po(x 2 / a) can be iden¬ 
tified as the gravitational potential due to the uniform Friedman background 
of density p 5 = po/a 3 . Hence the Lagrangian can be expressed as 
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is the difference between the total potential and the potential for the back 
ground Friedman universe &frw- Varying the Lagrangian in Eq.(3) with re¬ 
spect to x,;, we get the equation of motion to be: 
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Since 0 is the gravitational potential generated by the perturbed mass density, 
it satisfies the equation with the source (p — pb) = pbS: 

V^.0 = ATiGp b a 2 5 (6) 


Equation (5) and Eq.(6) govern the nonlinear gravitational clustering in an 
expanding background. 

Usually one is interested in the evolution of the density contrast 6 ( t , x) rather 
than in the trajectories. Since the density contrast can be expressed in terms 
of the trajectories of the particles, it should be possible to write down a dif¬ 
ferential equation for <5(f,x) based on the equations for the trajectories x(f) 
derived above. It is, however, somewhat easier to write down an equation for 
<5 k (t) which is the spatial Fourier transform of <5(t,x). To do this, we begin 
with the fact that the density p(x, t ) due to a set of point particles, each of 
mass m, is given by 

777 _ 

P(M) = H<M x - x *0)] ( 7 ) 


where Xj(f) is the trajectory of the ith particle and 5d is the Dirac delta 
function. The density contrast <5(x, t) is related to p(x, t) by 

1 + S(x, t) = x,;(t)] = J dq<Sr>[x-x T (t,q)]. (8) 


In arriving at the last equality we have taken the continuum limit by: (i) 
replacing x 2 ;(t) by xr(i, q) where q stands for a set of parameters (like the 
initial position, velocity etc.) of a particle; for simplicity, we shall take this to 
be initial position. The subscript l T’ is just to remind ourselves that xr(f, q) is 
the trajectory of the particle, (ii) replacing (V/N) by d 3 q since both represent 
volume per particle. Fourier transforming both sides we get 

<5 k (f) = J d 3 xe _lkx <5(x, t) = J d 3 q exp[—ik.x T (i, q)] — (27r) 3 (5 D (k) (9) 


Differentiating this expression, and using Eq. (5) for the trajectories give one 
can obtain an equation for <5 k (see e.g. ref. |9j; Eq.10). The structure of this 
equation can be simplified if we use the perturbed gravitational potential (in 
Fourier space) 0 k related to <5 k by 


fc 2 0k 

A'KGpba 2 


( k 2 a \ 
\47rGpoJ 


0k 



( 10 ) 


4 







In terms of 0k the exact evolution equation reads as: 


> k + 4-0 k = 
a 



zk.x 



(ii) 


where x = xr(t, q). Of course, this equation is not ‘closed’. It contains the 
velocities of the particles xy and their positions explicitly in the second term 
on the right and one cannot — in general — express them in simple form in 
terms of 0k- As a result, it might seem that we are in no better position than 
when we started. I will now motivate a strategy to tame this term in order to 
close this equation. This strategy depends on two features: 

(1) First, extremely nonlinear structures do not contribute to the right hand 
side of Eq.(ll) though, of course, they contribute individually to the two 
terms. More precisely, the right hand side of Eq.(ll) will lead to a density 
contrast that will scale as 5k oc k 2 if originally — in linear theory — 5^ oc k n 
with n > 2 as k —> 0. This leads to a P oc k 4 tail in the power spectrum fl, 1|. 
(We will provide a derivation of this result in the next section.) 

(ii) Second, we can use Zeldovich approximation to evaluate this term, once 
the above fact is realized. It is well known that, when the density contrasts 
are small, it grows as 5k oc a in the linear limit. One can easily show that such 
a growth corresponds to particle displacements of the form 

x T (a, q) = q - aV0(q); 0 = (47rGpo) _1 0 (12) 


A useful approximation to describe the quasi linear stages of clustering is 
obtained by using the trajectory in Eq.(12) as an ansatz valid even at quasi 
linear epochs. In this approximation, (called Zeldovich approximation), the 
velocities x can be expressed in terms of the initial gravitational potential. 

We now combine the two results mentioned above to obtain a closure condition 
for our dynamical equation. At any given moment of time we can divide the 
particles in the system into three different sets. First, there are particles which 
are already a part of virialized cluster in the non linear regime. Second set is 
made of particles which are completely unbound and are essentially contribut¬ 
ing to power spectrum at the linear scales. The third set is made of all particles 
which cannot be put into either of these two baskets. Of these three, we know 
that the first two sets of particles do not contribute significantly to the right 
hand side of Eq.(ll) so we will not incur any serious error in ignoring these 
particles in computing the right hand side. For the description of particles in 
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the third set, the Zeldovich approximation should be fairly good. In fact, we 
can do slightly better than the standard Zeldovich approximation. We note 
that in Eq. (12) the velocities were taken to be proportional to the gradient 
of the initial gravitational potential. We can improve on this ansatz by taking 
the velocities to be given by the gradient of the instantaneous gravitational 
potential which has the effect of incorporating the influence of particles in 
bound clusters on the rest of the particles to certain extent. 


Given this ansatz, it is straightforward to obtain a closed integro-differential 
equation for the gravitational potential. Direct calculation shows that the grav¬ 
itational potential is described by the closed integral equation: (The details of 
this derivation can be found in ref. [ 9 ] and will not be repeated here.) 


> k + 4-0 k = 


1 

3a 2 


_^P , , 

(27r) 3< ^ k+p< ^ k -P 




(13) 


This equation provides a powerful method for analyzing non linear clustering 
since estimating Eq.(ll) by Zeldovich approximation has a very large domain 
of applicability. In the next two sections, I will use this equation to study the 
transfer of power in gravitational clustering. 


3 Inverse cascade in non linear gravitational clustering: The k 4 
tail 


There is an interesting and curious result which is characteristic of gravita¬ 
tional clustering that can be obtained directly from our Eq. (13). Consider an 
initial power spectrum which has very little power at large scales; more pre¬ 
cisely, we shall assume that P{k ) oc k n with n > 4 for small k (i.e, the power 
dies faster than k 4 for small k). If these large spatial scales are described by 
linear theory — as one would have normally expected — then the power at 
these scales can only grow as P oc a 2 k n and it will always be sub dominant 
to k 4 . It turns out that this conclusion is incorrect. As the system evolves, 
small scale nonlinearities will develop in the system and — if the large scales 
have too little power intrinsically — then the long wavelength power will soon 
be dominated by the “fc 4 -tail” of the short wavelength power arising from the 
nonlinear clustering. This is a purely non linear effect which we shall now de¬ 
scribe. (This result is known in literature ||, @] but we will derive it from the 
formalism developed in the last section which adds fresh insight.) 

A formal way of obtaining the k 4 tail is to solve Eq. (13) for long wavelengths; 
i.e. near k = 0. Writing 0 k = 0 + 0® + .... where 0 k ^ = 0 k L) is the time 
independent gravitational potential in the linear theory and 0 k ; is the next 
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Fig. 1. The transfer of power to long wavelengths forming a /c 4 tail is illustrated 
using simulation results. Power is injected in the form of a narrow peak at L = 8. 
Note that the y —axis is (A/a) so that there will be no change of shape of the power 
spectrum under linear evolution with A oc a. As time goes on a k 4 tail is generated 
purely due to nonlinear coupling between the modes. (Figure adapted from ref. 00 

order correction, we get from Eq. (13), the equation 
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where £?(k, p) = [(7/8)fc 2 + (3/2)p 2 — 5(k-p/fc)) 2 . The solution to this equation 
is the sum of a solution to the homogeneous part [which decays as </> oc a -4 oc 
t~ 8 / 3 giving </> oc t ~ 5 / 3 ] and a particular solution which grows as a. Ignoring 
the decaying mode at late times and taking 0 k = aC k one can determine C k 
from the above equation. Plugging it back, we find the lowest order correction 
to be, 
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which is independent of k to the lowest order. Correspondingly the power 
spectrum for density Ps(k) oc a 2 k A P ip (k ) oc a A k A in this order. 

The generation of long wavelength k A tail is easily seen in simulations if one 
starts with a power spectrum that is sharply peaked in |k|. Figure 1 (adapted 
from 0) shows the results of such a simulation. The y-axis is [A (k)/a(t)\ 
where A 2 {k) = k 3 P/2n 2 is the power per logarithmic band in k. In linear 
theory A oc a and this quantity should not change. The curves labelled by 
a = 0.12 to a = 20.0 show the effects of nonlinear evolution , especially the 
development of k A tail. (Actually one can do better. The formation can also 
reproduce the sub-harmonic at L ~ 4 seen in Fig. 1 and other details; see ref. 

ft) 


4 Analogue of Kolmogorov spectrum for gravitational clustering 


If power is injected at some scale L into an ordinary viscous fluid, it cascades 
down to smaller scales because of the non linear coupling between different 
modes. The resulting power spectrum, for a wide range of scales, is well ap¬ 
proximated by the Kolmogorov spectrum which plays a key, useful, role in the 
study of fluid turbulence. It is possible to obtain the form of this spectrum from 
fairly simple and general considerations though the actual equations of fluid 
turbulence are intractably complicated. Let us now consider the corresponding 
question for non linear gravitational clustering. If power is injected at a given 
length scale very early on, how does the dynamical evolution transfers power 
to other scales at late times? In particular, does the non linear evolution lead 
to an analogue of Kolmogorov spectrum with some level of universality, in the 
case of gravitational interactions? 


Surprisingly, the answer is “yes”, even though normal fluids and collisionless 
self gravitating particles constitute very different physical systems. If power is 
injected at a given scale L = 2ir/ko then the gravitational clustering transfers 
the power to both larger and smaller spatial scales. At large spatial scales the 
power spectrum goes as P{k ) oc k A as soon as non linear coupling becomes 
important. We have already seen this result in the previous section. More in¬ 
terestingly, the cascading of power to smaller scales leads to a universal pattern 
at late times just as in the case of fluid turbulence. This is because, Eq.(13) 
admits solutions for the gravitational potential of the form cf^ft) = F(t)D( k) 
at late times when the initial condition is irrelevant; here F(t) satisfies a non 
linear differential equation and D(k) satisfies an integral equation. One can 
analyze the relevant equations analytically as well as verified the conclusions 
by numerical simulations. This study (the details of which can be found in 
ref. ll|) confirms that non linear gravitational clustering does lead to a uni¬ 
versal power spectrum at late times if the power is injected at a given scale 



initially. (In cosmology there is very little motivation to study the transfer of 
power by itself and most of the numerical simulations in the past concentrated 
on evolving broad band initial power spectrum. So this result was missed out.) 


Our aim is to look for late time scale free evolution of the system exploiting 
the fact Eq.(13) allows self similar solutions of the form </k(t) = F(t)D( k). 
Substituting this ansatz into Eq. (13) we obtain two separate equations for 
F(t) and -D(k). It is also convenient at this stage to use the expansion factor 
a(t ) = (t/fo) 2 ^ 3 of the matter dominated universe as the independent variable 
rather than the cosmic time t. Then simple algebra shows that the governing 
equations are 


d 2 F 7 dF 
da 2 ^ 2 da 


(17) 


and 


H 2 0 D k = 
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d 3 p 

(2^p 
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(18) 


Equation (17) governs the time evolution while Eq. (18) governs the shape of 
the power spectrum. (The separation ansatz, of course, has the scaling freedom 
F —> jiF^D —> (1 / g)D which will change the right hand side of Eq. (17) to 
— /aF 2 and the left hand side of Eq. (18) to /iFL^D-^. But, as to be expected, our 
results will be independent of yu; so we have set it to unity). Our interest lies 
in analyzing the solutions of Eq. (17) subject to the initial conditions -F(aj) = 
constant, ( dF/da)i = 0 at some small enough a = at. 


Inspection shows that Eq. (17) has the exact solution F(a) = (3/2)a _1 . This, of 
course, is a special solution and will not satisfy the relevant initial conditions. 
However, Eq. (17) fortunately belongs to a class of non linear equations which 
can be mapped to a homologous system. In such cases, the special power 
law solutions will arise as the asymptotic limit. (The example well known 
to astronomers is that of isothermal sphere □ Our analysis below has a 
close parallel.) To End the general behaviour of the solutions to Eq. (17), we 
will make the substitution F(a) = (3/2 )a~ 1 g(a) and change the independent 
variable from a to q = log a. Then Eq. (17) reduces to the form 


d 2 g 1 dg 3 

dy + 2^ + 2 s(j - 1) = 0 


(19) 


This represents a particle moving in a potential V(g) = (1/2 )g 3 — (3/4 )g 2 
under friction. For our initial conditions the motion will lead the “particle” 
to asymptotically come to rest at the stable minimum at g — 1 with damped 
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Fig. 2. The solution to Eq. (19) plotted in g — a plane. The function g(a) asymp¬ 
totically approaches unity with oscillations which are represented by the spiral in 
the right panel. The different curves in the left panel corresponds to the rescaling 
freedom in the initial conditions. One fiducial curve which was used to model the 
simulation is shown in the red. For more details, see ref. Jll| . 

oscillations. In other words, F(a) —> (3/2)a _1 for large a showing this is indeed 
the asymptotic solution. From the Poisson equation, it follows that k 2 (p k oc 
(hk/o) so that 5k(a) oc g(a)k 2 D( k) giving a direct physical meaning to the 
function g(a) as the growth factor for the density contrast. The asymptotic 
limit (g ~ 1) corresponds to to a rather trivial case of 5k becoming independent 
of time. What will be more interesting — and accessible in simulations — will 
be the approach to this asymptotic solution. To obtain this, we introduce the 
variable 


u = g + 2 



( 20 ) 


so that our system becomes homologous. It can be easily shown that we now 
get the first order form of the autonomous system to be 

du = 6 g(g - 1) ^ 

dg u-g 


The critical points of the system are at (0,0) and (1,1). Standard analysis 
shows that: (i) the (0, 0) is an unstable critical point and the second one (1,1) 
is the stable critical point; (ii) for our initial conditions the solution spirals 
around the stable critical point. 


10 







Fig. 3. Left panel: The results of the numerical simulation with an initial power 
spectrum which is a Gaussian peaked at L = 24. The y-axis gives A& where 


A? = k?P/2'K 2 is the power per logarithmic band. The evolution generates a well 


known k 4 tail at large scales (see for example, jlflj) and leads to cascading of power 
to small scales. Right panel: The simulation data is re-expressed by factoring out 
the time evolution function g(a ) obtained by integrating Eq. (19). The fact that the 
curves fall nearly on top of each other shows that the late time evolution is scale free 
and described by the ansatz discussed in the text. The rescaled spectrum is very 
well described by P(k) oc k~ 0A /( 1 + (fc/fco) 2 ' 6 ) which is shown by the, completely 


ul- 


overlapping, broken blue curve. For more details, see ref. 

Figure 2 (from ref. Tlj) describes the solution in the g — a plane. The g(a) 
curves clearly approach the asymptotic value of g cs 1 with superposed oscil¬ 
lations. The different curves in Fig. (2) are for different initial values which 
arise from the scaling freedom mentioned earlier. (The thick red line corre¬ 
spond to the initial conditions used in the simulations described below.). The 
solution g(a) describes the time evolution and solves the problem of determin¬ 
ing asymptotic time evolution. 


To test the correctness of these conclusions, we performed a high resolution 


simulation using the TreePM method |13|, [1J] and its parallel version 


with 128 3 particles on a 128 3 grid. Details about the code parameters can be 


found in 14]. The initial power spectrum P(k) was chosen to be a Gaussian 


peaked at the scale of k p = 2i t/L p with L p = 24 grid lengths and with a 
standard deviation A k = 2n/Lbox, where Lb ox = 128 is the size of one side 
of the simulation volume. The amplitude of the peak was taken such that 
A Hn (k p = 27t/ L p , a = 0.25) = 1. 


The late time evolution of the power spectrum (in terms of A 2 . = k 3 P(k)/2n 2 
where P = |<5fc| 2 is the power spectrum of density fluctuations) obtained from 
the simulations is shown in Fig.3 (left panel). In the right panel, we have 
rescaled the A fc , using the appropriate solution g(a ). The fact that the curves 
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fall on top of each other shows that the late time evolution indeed sales as 
g(a ) within numerical accuracy. A reasonably accurate fit for g(a ) at late 
times used in this figure is given by g(a) oc a(l — 0.3 In a). The key point to 
note is that the asymptotic time evolution is essentially 5(a) oc a except for a 
logarithmic correction, even in highly nonlinear scales. (This was first noticed 
from somewhat lower resolution simulations in H3o- Since the evolution at 
linear scales is always 5 oc a, this allows for a form invariant evolution of 
power spectrum at all scales. Gravitational clustering evolves towards this 
asymptotic state. 

To the lowest order of accuracy, the power spectrum at this range of scales is 
approximated by the mean index n ~ — 1 with P(k) oc /c -1 . A better fit to the 
power spectrum in Fig.3 is given by 

k~ 0A 2vr 

l + (k/ko)™' h 




This fit is shown by the broken blue line in the figure which completely overlaps 
with the data and is barely visible. (Note that this fit is applicable only at 
L < L p since the fc 4 tail will dominate scales to the right of the initial peak; 
see the discussion in H3). At nonlinear scales P(k) oc k~ 3 making A k flat, as 
seen in Fig. 3. (This is not a numerical artifact and we have sufficient dynamic 
range in the simulation to ascertain this.) At quasi linear scales P(k) oc k~ 0A . 
The effective index of the power spectrum varies between —3 and —0.4 in this 
range of scales. 


5 Possible interpretation 


In the case of viscous fluids, the energy is dissipated at the smallest scales as 
heat. In steady state, energy cannot accumulate at any intermediate scale and 
hence the rate of flow of energy from one scale to the next (lower) scale must 
be a constant. This constancy immediately leads to Kolmogorov spectrum. 
In the case of gravitating particles, there is no dissipation and each scale 
will evolve towards virial equilibrium. At any given time t, the power would 
have cascaded down only up to some scale l m in(t ) which it self, of course, is 
a decreasing function of time. So, at time t we expect very little power for 
1 < kl min (t ) and a k A tail for kL p < 1, say. The really interesting band is 
between l min and L p . 


To understand this band, let us recall that the Lagrangian in Eq.(3) leads 
to the time dependent Hamiltonian is H(p, x, t) = J2[p 2 /2rna 2 + U}. The 
evolution of the energy in the system is governed by the equation dH/da = 
(dH/da) PiX . It is clear from Eq. (4) that (dU/da) PtX = —[//awhile (<9T/<9a) PiX = 
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—2 T/a. Hence the time evolution of the total energy H = E of the system is 
described by 


d.E 1 m rrN 2 E U E T 

— = — 2 T + U) = -+ - =- 

da a a a a a 


(23) 


In the continuum limit, ignoring the infinite self-energy term, the potential 
energy can be written as: 


= _M f d?* f'p y S( -*; a)S(y ; a) 

2a J J x — y 


(24) 


Hence 

&- = aU = -M /> x 

da 2 J J x — y 


(25) 


The ensemble average of U, per unit volume, will be 


£ 


Gpt 

2 Va 


dW y {Hx ; a)Siy ; a)) 

l x — y| 


oc 



(26) 


where V is the comoving volume. 

When a particular scale is virialized, we expect £ ~ constant at that scale in 
comoving coordinates. That is, we would expect the energy density in Eq.(26) 
would have reached equipartition and contribute same amount per logarithmic 
band of scales in the intermediate scales between l min and L pea k. This requires 
P(k) oc 1 /k which is essentially what we found from simulations. The time 
dependence of P is essentially P oc a 2 (except for a logarithmic correction) in 
order to maintain form invariance of the spectrum with respect to the linear 
end. Similarly the scale dependence is P oc k~ l which is indeed a good fit to 
the simulation results. The flattening of the power at small scales, modeled by 
the more precise fitting function in Eq.(22), can be understood from the fact 
that, equipartition is not yet achieved at smaller scales. The same result holds 
for kinetic energy if the motion is dominated by scale invariant radial flows 
BE3- Our result suggests that gravitational power transfer evolves towards 
this equipartition. 
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